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Abstract. An algorithm for linear estimation of aerosol bulk 
properties such as particle volume, effective radius and com- 
plex refractive index from multiwavelength lidar measure- 
ments is presented. The approach uses the fact that the total 
aerosol concentration can well be approximated as a linear 
combination of aerosol characteristics measured by multi- 
wavelength lidar. Therefore, the aerosol concentration can 
be estimated from lidar measurements without the need to 
derive the size distribution, which entails more sophisticated 
procedures. The definition of the coefficients required for 
the linear estimates is based on an expansion of the particle 
size distribution in terms of the measurement kernels. Once 
the coefficients are established, the approach permits fast re- 
trieval of aerosol bulk properties when compared with the 
full regularization technique. In addition, the straightforward 
estimation of bulk properties stabilizes the inversion making 
it more resistant to noise in the optical data. 

Numerical tests demonstrate that for data sets containing 
three aerosol backscattering and two extinction coefficients 
(so called 3/1 + 2 a) the uncertainties in the retrieval of parti- 
cle volume and surface area are below 45 % when input data 
random uncertainties are below 20 %. Moreover, using linear 
estimates allows reliable retrievals even when the number of 
input data is reduced. To evaluate the approach, the results 
obtained using this technique are compared with those based 
on the previously developed full inversion scheme that re- 
lies on the regularization procedure. Both techniques were 
applied to the data measured by multiwavelength lidar at 
NASA/GSFC. The results obtained with both methods using 


the same observations are in good agreement. At the same 
time, the high speed of the retrieval using linear estimates 
makes the method preferable for generating aerosol infor- 
mation from extended lidar observations. To demonstrate the 
efficiency of the method, an extended time series of obser- 
vations acquired in Turkey in May 2010 was processed us- 
ing the linear estimates technique permitting, for what we 
believe to be the first time, temporal-height distributions of 
particle parameters. 


1 Introduction 

Theoretical and experimental studies of the last decade 
have demonstrated that multiwavelength (MW) Raman li- 
dars based on a tripled Nd:YAG laser are able to provide 
the height distribution of particle physical parameters, such 

as radius, concentration and complex refractive index (Ans- 
mann and Muller, 2005). Moreover, up to a certain limit 

such systems can reproduce the main features of the parti- 
cle size distribution in the 0.075-10 pm radii range. To in- 
vert the aerosol extinction a and backscattering /J coeffi- 
cients measured at multiple wavelengths to particle param- 
eters, numerous possible approaches have been considered 
but for routine processing of lidar measurements the inver- 
sion with regularization is now the most commonly used (see 
Muller et al., 1999; Veselovskii et ah, 2002, 2004, 2009; Kol- 
gotin and Muller, 2008, and references therein). In order to 
adequately address the fundamental non-uniqueness of the 
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lidar data interpretation, a family of solutions is generated in 
the framework of this approach. Specifically, a series of so- 
lutions is generated using different initial guesses, different 
aerosol assumptions and different settings of a priori con- 
straints. Each single solution is obtained using the regulariza- 
tion technique. Then the individual solutions corresponding 
to the smallest residuals are averaged and the result of the av- 
eraging is taken as the best estimate of the aerosol properties. 

This approach has demonstrated the possibility of provid- 
ing rather adequate retrievals of aerosol properties. How- 
ever it is quite time-consuming, a fact that becomes an is- 
sue when large volumes of data need to be analyzed, as for 
example from an air- or space-borne lidar system. Installa- 
tion of MW lidars on air or space-borne platforms poses an- 
other problem: the retrieval algorithm should be more tol- 
erant to noise in the input data since reasonable averaging 
times are likely to be smaller for moving lidar systems. And 
finally, in the regularization approach described in (Muller 
et al., 1999; Veselovskii et al., 2002) at least five input opti- 
cal data (three backscatterings and two extinctions, so called 
3 P + 2a) are needed to retrieve the particle size distribution 
(PSD), but in many applications it would be highly desir- 
able to decrease the number of optical channels. So, the 
development of the approach permitting the reliable esti- 
mation of particle bulk properties such as volume, surface 
density and effective radius from a reduced number of opti- 
cal channels would be an important improvement. One way 
to assess this possibility is to attempt to approximate the 
bulk properties by a linear combination of the input opti- 
cal data (extinction and backscattering). The correspond- 
ing weight coefficients can be determined by expanding the 
PSD in terms of the measurement kernels (Twomey, 1977). 
Thomason and Osborn (1992) used this approach to estimate 
aerosol mass with the multiwavelength SAGE II extinction 
kernels. The interpretation of lidar measurements using the 
linear estimate techniques was explored in early studies by 
Chaikovskii and Shcherbakov (1985). The potential of this 
approach for treating elastic-Raman multiwavelength lidar 
measurements was studied by Donovan and Carswell (1997) 
under assumption of known refractive index. The technique 
was further explored in recent publications (De Graaf et al., 
2009, 2010) where different aerosol models were used to in- 
vert optical data without prior information about the particle 
refractive index. 

In this paper we propose a modified technique, which 
here and below we refer to as “linear estimation" (LE). In 
difference with the mentioned above approaches the com- 
plex refractive index is derived as a part of retrieval proce- 
dure together with the bulk aerosol characteristics. In ad- 
dition, in order to improve stability of solution we provide 
not a single solution but a family of solutions closely repro- 
ducing the measurements. To validate LE, we apply it and 
the full inversion algorithm (Veselovskii et al., 2009) to the 
same data and compare the results. Finally, we apply LE to 


an extended series of lidar measurements to evaluate height- 
temporal variations of the particle bulk parameters. 

2 Algorithm description 

The aerosol extinction (a) and backscattering coefficients 
(fi) are related to the particle volume size distribution v ( r ) 
via integral equations as follows: 

f max 

Sp= f Ki(m,r) v(r) dr l = (i, X k ) = 1, . . ., No (1) 

r min 

Index l labels the type of optical data (i — a, fi) and wave- 
lengths A. jt; Kfim , r) are the volume kernels (VK) depending 
on the complex refractive index m — ;«r — i ■ m\ and particle 
radius r e[r m i n ,r max ]. To get kernels Ki(m,r) in our study 
the Mie computations are used, thus the particles are as- 
sumed to be spherical. This approach can also be generalized 
to treat the particles of irregular shape, by using the kernels 
corresponding the ensemble of randomly oriented spheroids 
(Mishchenko et al., 2000; Dubovik et al., 2006). In vector- 
matrix form Eq. (1) can be rewritten as: 

S = Ku (2) 

Here v is the column vector with elements Vk corresponding 
to the particle volume inside radii interval [r>, 1 ] and K 

is the matrix containing the discretized kernels as rows. The 
volume distribution v(r) can be expanded in terms of the ker- 
nels of Eq. (1), as prescribed by Twomey (1977). Such an ex- 
pansion assumes that vector v corresponding to v(r) can be 
presented as a combination of the matrix K rows, i.e. 

v = v g + v± — K t jc + Dj_ (3) 

where v g is the projection of the volume distribution on the 
measurement kernels, while v± is the residual - the part 
of volume distribution orthogonal to these kernels (Kuj_ = 
0) and xj are the weight coefficients of expansion. Using 


Eq. (3), Eq. (2) can be rewritten as: 

g — KK t jc + K i>i = KK t x (4) 

Then, the vector x of the expansion coefficients can be found 
as: 

x=(KK T )“T (5) 

and the volume distribution projected on the kernels is 
v g = K t x = K t (KK t ) ” 1 g (6) 

The residual then can be written as: 

I)_L = V - Vg = (i - K t (kK t ) 1 Kj v (7) 
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Equations (6)-(7) can now be used to evaluate the linear esti- 
mations of the unmeasured aerosol characteristics. If there is 
one or several aerosol characteristics p, ( i — 1, . . . , N p ) that 
are not measured but are needed to be estimated using mea- 
surements g , the dependence of /?, on the size distribution 
can be described as: 

p = Pv. (8) 

Here the elements p, of vector p are the unmeasured aerosol 
characteristics, and P is the matrix of the corresponding co- 
efficients. Taking into account Eq. (3) the vector p can be 
expressed as: 

p = P(v g + u_l) = P g + Pj_=^g + Dj_i> (9) 

Here p„ represents the vector of projections of characteris- 
tics pi on the measured set g and p represents the vector 
of characteristics pt on the null-space t>j_. In other words, p„ 
can be estimated from g , while the measurements provide no 
information about p . Using Eqs. (6) and (7) the matrices of 
coefficients F and 1)^ can be expressed as 

F =PK T (KK T )~ 1 and 

Dj_ = P ^7 - K t (KK t ) 1 (10) 

The elements of matrix F can be computed and stored in 
the look-up tables making computations of p g very fast. The 
residual term p j_ cannot be measured with the available set of 
observations g, but can be estimated from numerical model- 
ing for typical situations. The situation is particularly favor- 
able when particle bulk property p (for example, volume, 
surface, number density) needs to be estimated (Donovan 
and Carswell, 1997). In this case the matrix P contains the 
weight coefficients for different integral properties as rows. 
For example, for volume (i = 1) P\t — 1, for surface (i — 2) 
P 2 k = y: and for number density (i — 3) P}k = ■ In such 

case, the existence of the zero space v± does not have much 
importance and the residual p is generally expected to be 
small, because the observations g are known to be strongly 
sensitive to aerosol total concentrations, while being less 
sensitive to the details of the size distribution. 

Thus, the projection p g can be estimated quickly from the 
observations g without calculating the full size distribution 
v, i.e. without performing a full inversion of Eq. (2). This 
is a significant advantage of the using the linear estimates 
p, ; as compared to the more conventional approach which 
yields PSD (e.g. Muller et al., 1999; Veselovskii et al., 2002, 
2004, 2009; Kolgotin and Muller, 2008). Indeed, the calcu- 
lation of p q is fast and defining the coefficients F is straight- 
forward. Following Eq. (10) the calculation of F involves the 
inversion of matrix KK T . In principle this operation can be 
ambiguous if KK T is ill-conditioned. However, in particular 
cases when only a very few measured characteristics gj are 


used (for example, in our case we use maximum 5 different 
observations), the matrix KK T has small dimension (maxi- 
mum 5 x 5) and in the case such as here that each of the five 
measured characteristics g; is quite different, is well-posed 
and can be inverted exactly. By contrast, the conventional ap- 
proaches that provide the size distribution v and must solve 
Eq. (2) may face significant difficulties. For example, the 
least square solution can of Eq. (2) is: 

u=(k t k)“‘ K T g (11) 

Here the matrix K T K has to be inverted. This matrix has di- 
mension N v x N v , where N v is the dimension of size dis- 
tribution v. Therefore, matrix K T K has significantly larger 
dimension than matrix KK T . For example, the aerosol 
retrievals from sun-photometer observations discussed by 
Dubovik and King (2000) use N v = 22. In such situations 
K t K is known to be ill-conditioned and the inversion of 
this matrix becomes ambiguous. Therefore in many practi- 
cal applications different types of constraints can be used to 
achieve unique and stable solution of Eq. (2). For example, 
it can be constrained based on smoothness of the solution as 
suggested by Phillips (1962) and Twomey (1977). However, 
the use of smoothness or other a priori constraints may re- 
quire rather sophisticated developments (see discussion by 
Dubovik, 2004). We should mention also that in most of al- 
gorithms for lidar data inversion (e.g. Muller et al., 1999; 
Veselovskii et al., 2002, 2004, 2009; Kolgotin and Muller, 
2008), the size distribution is described by a much smaller 
number of parameters than 22 (typically N v = 5-7). This re- 
duces the difficulties associated with the inversion of matrix 
K t K, however this also may lead to the introduction of ad- 
ditional errors in the algorithm, since some features of the 
size distribution are neglected. In this respect the coefficients 
F can be calculated using the detailed size distribution with 
very large N v , since matrix KK T has the same dimension as 
number of rows of K. Correspondingly, coefficients F can 
be always found accurately. In our computations we nor- 
mally use N v = 100 radii logarithmically distributed inside 
the inversion interval. 

The measured optical data g* contain the error A g 
g* = g + A ff (12) 

Thus the uncertainty of particle parameters estimation is: 

Ap = Fg* - p = F (g + A g ) - (F g + D ± i;) 

= FA g -D ± i> (13) 

Then, if the measurement errors A g are random, unbi- 
ased (i.e. (Ag) = 0) and have covariance matrix I^A g A^j — 
C„, the corresponding covariance matrix of retrieval 

uncertainties A p can be written as 

C p = ((FA* - Dj_d) (FA g - Dj_d) T | 

£ random _|_ ^.systematic ( 14 ) 
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Here, term C™ 1110 " 1 represents the contribution of the random 

measurement errors A g to C p and c^ ystematic represents the 
non-random part of the errors appearing due to existence of 
»i which is orthogonal to the kernels K. These terms can be 
expressed as: 

c random = F | Ag A^F T = FC ? F T 

= PK T (KK T ) ' C 4 , (kk'J KP t (15) 

c systematic = p±p T = D± 

= P - K t (KK T ) 1 K ^ (vv T ^j 

-K t (kK t ^ _ 1 K^P t (16) 

The first term can be calculated using the covariance matrix 
of the measurements CL, and the second term can be esti- 
mated using a priori estimates of v. For example, our model- 
ing experiments in the next section show that p have gen- 
erally rather small values and, therefore, the accuracy of the 
retrieval of p performed using the “projection” p g (see Eq. 9) 
is acceptable in the majority of cases. 

All equations given above are written with the assump- 
tion that the matrix K in Eq. (2) and, therefore, matrix F 
in Eq. (10) are known accurately. However, this is not really 
the case since K depends on the complex refractive index. 
Therefore, if the actual value of the complex refractive index 
m(X) is not known and we use an estimate m(X), instead of 
Eq. (10) we have: 

p = P(v g + v ± ) = PK T (KK T y l g + Pi ± (17) 

or it can be written in the same manner as Eq. (9) 

P = Pg + P± = Fg + D±» (18) 

Correspondingly, if we use an estimate of the complex re- 
fractive index, the estimate p* = F g*, should be replaced by 

p* = Fg*. The corresponding uncertainties of the retrieval 
can be estimated as: 

A P = P* g ~P 

= Fg* - (P g + Pa) = F(g + A g ) - (Fg + Dj_i>) 

= FAg — AFg — Dj_u = FAg — AFKu — Dj_c 
= FAg — (AFK + DjJu (19) 

where AF = F — F. The covariance matrix of uncertainties 
is: 

C p = ((FAg-(AFK + Dj>) 

(FAg-(AFK + D ± )t>) T 

^random , ^systematic 

— \^ p -r 


where 

c random = f A = FCgF T 

= PK T (M T j ' Cg (kK T ) 'kP T (21) 


c systematic = ( AFK + D±) ( A FK + D ± ) T 

= AFK K t (AF) t + AFK 

+D ± ( v t> T ) K t ( AF) t + D ± ( v u T ) (22) 


Thus, if we compare Eqs. (15) and (21), the only difference 
is that Eq. (21) uses matrix of coefficients F instead of F. 
It is reasonable to expect that the magnitudes of elements 
of F are close to those of F since the optical characteristics 
generally do not exhibit very high sensitivity to variations 
of m(X). Therefore one can expect that the components of 
£ random gj ven by gq (21) will have magnitudes close to those 
given by Eq. (15). 

By contrast C p in Eq. (22) comparing to Lf, 

of Eq. (16) has three extra terms containing AF. If AF = 0 
Eq. (22) coincides with Eq. (16). Another observation is that 
if we have a rather complete set of observations g*, so that we 
do not have a null-space, i.e. Ily = 0, then Eq. (14) retains 
only one first term C™ Klom , while Eq. (20) still retains the 
second term that represents the systematic bias: 


C „ = C p 


random 


^.systematic 


AFK (™ T ) K T ( AF) 


— FC ? F 

T 


(23) 


Thus, the use of estimate m(X) that is different from the ac- 
tual value of complex refractive index m (A.) always leads to 
the appearance of a systematic error term C^) stemalic . More- 
over, it is rather clear that c^ ystemauc may easily dominate 
over C' ; ;‘ ndom in Eq. (23). This becomes rather obvious when 
using g = Ku in Eq. (23): 


C p = FCgF T + AF(gg T )(AF) T 


(24) 


Here the first term contains C g = [A g A g j - the covariance 
matrix of errors A g of g and the second term contains the 
matrix gg T . Since the magnitude of errors A g are generally 
much smaller than the magnitudes of measurements g, the 
elements of Cg are much smaller than gg T . Therefore, even 
if AF = F — F is not very significant the magnitude of the 
second term of Eq. (23) is likely to remain considerable. 

Thus, if the sensitivity of g* to m ( X ) is high, the errors due 
to wrongly chosen m(X) can be much higher than errors due 
to measurement uncertainties and existence of null-space. At 
the same time if the sensitivity of g* to m(X) is high, and set 
of observations (g*) T = ^gjf; gp , ...; g/y 0 ) is quite represen- 
tative we can attempt to estimate m(X) from available obser- 
vations. The input optical data (backscatters and extinctions) 
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are themselves the particle properties and each g* can be re- 
calculated back from the rest of Ng- 1 data using Eq. (9), as 
suggested in (De Graaf et al., 2009, 2010). By doing so for 
each optical data, we get Ng estimates of gj that we compare 
with the observations g*j. It should be mentioned that we can 
not make these estimates using all Ng optical data, because in 
that case gj and g* coincide. If the sensitivity of g* to m (X) is 
high the magnitudes of errors A gj — gj — g*. should strongly 
depend on the assumed value of m(X). Correspondingly, if 
we obtained the set of estimates gj (m) using different as- 
sumed values of m(X) we can attempt to estimate m(X) by 
searching for the smallest errors A gj — gj — g for exam- 
ple by searching for the minimum of the following quadratic 
form: 

Vk (7 n) = (g (777) - g*) T C” 1 (Ir (m) - g*) (25) 

If the sensitivity of g* to m (X) is high this form should have a 
well-defined minimum and 777(A) can be estimated using the 
available measurements g*. In our study we assume that the 
errors of the measurements are the same for all channels, and 
refractive index is spectrally independent inside the spectral 
range that is considered. Then the refractive index is found 
from the minimum of discrepancy: 

e (*;-«(»)) 2 


Since there is no a priori knowledge about the particle size 
distribution and refractive index we find the discrepancy p 
for all predefined values of 7' m i n , r max , lying in the inter- 
val 0.075 pm-10 pm, and for the set of values m r and mi 
from respective intervals 1.35-1.65 and 0.00-0.03 just as we 
did in our regularization algorithm (Veselovskii et al., 2002). 
Normally the total number of predefined combinations does 
not exceed Nj = 3000. Based on our previous experience 
(Veselovskii et al., 2002) we prefer to average the solutions 
near the minimum of discrepancy rather than take a single so- 
lution. Such an averaging procedure stabilizes the inversion. 
To choose the averaging interval the solutions are ranged in 
accordance with their discrepancy from minimum p mm to 
maximum p max , normally 1 % of solutions are averaged. The 
retrieval for each vertical bin was done independently. 

The estimation of the refractive index from the minimiza- 
tion of p in Eq. (26) is illustrated by Fig. 1. The discrepancy 
p and the uncertainty of the volume retrieval ey are given 
for different assumptions of the value of 7?7r. Synthetic input 
data were generated assuming a log-normal aerosol distribu- 
tion with modal radius 7 - 0 = 1 pm and variance 0.4; the 
model refractive index is m — 1.5-/0.005. Both p and ey 
have minima at 777 R =1.5 corresponding to the “true” value 
of 777 r, thus the minimization of the discrepancy minimizes 
the uncertainty of the particle volume estimation. Similar 
plots can be provided also for the imaginary part of the re- 
fractive index. The presence of noise A g in the input data 



Fig. 1. Dependence of discrepancy p and uncertainty of particle vol- 
ume estimation sy on the real part of refractive index. Simulation 
was performed for rg = 1 pm and m = 1.5-/0.005. 

naturally increases the minimum value of discrepancy p that 
can be achieved and thus the uncertainty of the refractive in- 
dex retrieval. The actual increase of the retrieval uncertainly 
also depends on the particle size distribution and specific re- 
alization of errors A g in the optical data. To evaluate the cor- 
responding uncertainties of the estimation of particle param- 
eters numerical simulations using different types of PSD and 
different input errors can be performed as we will illustrate 
in the next section. 

Thus the main difference of described in this section algo- 
rithm from the approach presented previously by Donovan 
and Carswell (1997); De Graaf et al. (2009, 2010), is that we 
consider not a single solution but a family of linear solutions 
corresponding different inversion intervals 7- m ; n , r max and dif- 
ferent complex refractive indices. The average of solutions 
in the vicinity of the minimum of discrepancy (Eq. 26) is 
considered as most probable estimate of particle parameters. 

3 Estimation of retrieval uncertainties 

Numerical simulation is used here to test the algorithm and 
to estimate the retrieval uncertainties. In these simulations 
we used synthetic input optical data assuming a log-normal 
aerosol distribution with rg — 0.2, and 2 pm, which 

are typical values for the fine and the coarse mode particles 
(e.g. see Dubovik et al., 2002), and the variance in all cases 
is 0.4. As discussed in the previous section, one of the prin- 
cipal questions in the application of the LE technique is the 
estimation of the residual v± in Eq. (3). This residual will de- 
pend on the PSD and the refractive index. The computations 
performed demonstrate that for all values of m considered 
here, the residual v± is below 4 % and 15 % for r g — 0.2, and 
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Fig. 2. Cumulative probability of uncertainty sy of volume density 
retrieval from 3/1 + 2 a data with input errors e = 10 %, 20%, 30 %. 
Simulation was performed for rg = 0.2 pm using volume kernels. 

2 pm respectively. So the existence of a null-space does not 
present a serious limitation to the LE technique for typical 
atmospheric aerosols. 

To evaluate the effect of input uncertainties, the random er- 
rors in the range of [ 0 , ±e] were added to the data and from 
these distorted optical data, the particle parameters were re- 
trieved. We assume that the uncertainties in all measurement 
channels are equivalent so that all the diagonal elements of 
the error covariance matrix C g = ^A g A^j are the same. The 
procedure was repeated 1000 times allowing robust statistics 
to be gathered. The retrieval uncertainties are presented in 
the form of probability distributions such as shown in Fig. 2 
where a typical cumulative probability of volume density un- 
certainty sy is shown. For every value of Ey the plot gives the 
probability that the retrieval uncertainty is below this value. 
For example, from the plots in Fig. 2 we can conclude that 
in 90 % of the cases the spread in the values of the volume 
estimation is below sy ~ 20 %, 35 %, 50 % for input errors 
e=10%, 20%, 30% respectively. We take these values to 
represent the uncertainty in the retrieval. Thus the uncer- 
tainty rises approximately linearly with £ and the method can 
provide reasonable estimations even for 30% input errors. 

The results shown in Fig. 2 were obtained using volume 
kernels (VK) of Eq. (1), but Eq. (1) can be written also us- 
ing other types of kernels corresponding to the number 
surface ^ or volume density size distribution in log- 
arithmic space. All these kernels (henceforth referred to as 
NK, SK and VLK ) can also be used in retrievals. Donovan 
and Carswell (1997), reported that in their approach for the 
retrieval of surface density the surface kernels were prefer- 
able, while for volume retrieval the volume kernels were 
better suited. In our study, we also tested different types of 
kernels. We did not notice a significant difference between 
these kernels for small particles, but for particles in the coarse 


Table 1. Uncertainties of particle parameters estimation. Results are 
obtained for PSDs with modal radii 0.2 pm and 2 pm for input 
errors e = 0, 10 %, 20 %. 


r 0> 1™ 

Input random 
uncertainties 


0.2 



2 


0 

10% 

20% 

0 

10% 

20% 

Sy, % 

5 

20 

35 

15 

30 

45 

e s , % 

5 

20 

45 

2 

10 

30 

% 

5 

25 

40 

15 

25 

35 

s N , % 

10 

40 

60 

25 

75 

110 

Sm R 

0.01 

0.05 

0.07 

0.015 

0.025 

0.04 


mode the volume kernels provided slightly better estimations 
of all parameters. The difference with the results of (Dono- 
van and Carswell, 1997) may be due to the optimization of 
inversion intervals in our algorithm. All retrievals presented 
below were obtained with the volume kernels. 

The simulation results are summarized in Table 1 showing 
the uncertainty of volume (£y), surface (£ 5 ), number ( sn ) 
density, effective radius (£« eff ), and real part of refractive in- 
dex (£ mR ) retrieval (taken at 90 % probability level) for input 
random uncertainties of £ = 0, 10%, 20%. The effective ra- 
dius was estimated from the ratio of volume and surface den- 
sity: r e ff = 3 j . The results are given for ro = 0.2 pm and 2 pm 
to separately characterize uncertainties for small and big par- 
ticles. In the absence of input errors, the uncertainties of the 
retrieval are due to the null-space and the unknown value of 
the refractive index as follows from Eq. (22). Minimization 
of discrepancy (Eq. 26) keeps the uncertainty of the volume 
estimation below 5 % for small particle sizes characteristic of 
the fine mode and below 15 % for particles with sizes more 
consistent with the coarse mode particles. 

From the results shown in Table 1, several conclusions can 
be made. First of all, the retrieval is stable for both small and 
big particles and an uncertainty of volume estimation below 
45 % can be obtained even for 20% input errors. For small 
particles the uncertainties of surface, volume and effective 
radius estimation are close, while for big particles the sur- 
face density is the most stable parameter in retrieval where 
the corresponding uncertainty of £5 is less than 30 % even 
for 20% input errors. The most unstable parameter in the 
retrieval is the number density where the corresponding un- 
certainty for particles with ro = 2 pm is above 100 % for 20 % 
input errors. The real part of particle refractive index can be 
retrieved more accurately for big particles, where the cor- 
responding uncertainty of the real part is below ±0.04 for 
£ = 20 %, while for small particles this uncertainty increases 
to ±0.07. 

In our retrievals we considered the full data set 3/1 + 2 a 
and the reduced one 3/1 + la, where extinction at 532 nm 
was removed. The important finding is that the uncertain- 
ties in the estimates of particle parameters from 3/1+1 a 
data in most cases did not exceed the corresponding values 
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for 3/3 + 2 a, thus it seems apparent that the number of in- 
put data can be decreased when only particle bulk proper- 
ties are desired. Evaluation of retrieval uncertainties for dif- 
ferent combinations of the optical data and different parti- 
cles characteristics is in our plans but beyond the scope of 
present paper. 

The imaginary part of the refractive index is one of the 
most difficult parameters to estimate from multi-wavelength 
lidar as the kernels are not very sensitive to changes in the 
value of m\. As already mentioned, the inverse problem (1) is 
strongly underdetermined, so the solution depends on the 
constraints used, in particular on the range of refractive index 
values considered during the minimization of the discrepancy 
(Eq. 26). To evaluate the influence of the range of m i consid- 
ered on the retrievals, we performed simulations for three in- 
tervals: 0 < m\ < 0.01, 0 < /«i < 0.02 and 0 < mi < 0.03 as- 
suming 10% errors in input data and model value m = 1.5- 
/0.005. Computations performed for particles with kq = 2 pm 
show that the uncertainties in the estimate of m\ for these 
intervals are 50%, 100% and 140% respectively. Hence, 
for reasonable estimation of the imaginary part of the re- 
fractive index it is very desirable to have a priori informa- 
tion about the aerosol type to constrain the range of mi that 
is considered. 

To illustrate the influence of ni\ on the estimation of other 
parameters. Fig. 3 shows the cumulative probability plots for 
the volume retrieval of particles with r {) = 0.2 pm and 2 pm 
using the three ranges of mi mentioned above. For particles 
in the coarse mode the uncertainty in sy increases from 25 % 
to 35 % when maximal value of mi rises from 0.01 to 0.03. 
For particles in the fine mode the retrieval is essentially in- 
sensitive to the range of m\ considered. Thus, in spite of the 
ambiguity in the retrieval of the imaginary part, the uncer- 
tainty in mi has little influence on the estimation of the other 
parameters. 

4 Comparison with regularization retrievals 

To validate the approach described in Sect. 2, the linear es- 
timation (LE) and regularization (Veselovskii et al., 2002) 
algorithms were applied to the same experimental data ob- 
tained by multiwavelength Raman lidar at NASA/GSFC in 
Greenbelt, MD during August-September 2006 (Veselovskii 
et al., 2009). The lidar is based on a tripled Nd:YAG laser 
and provided three particle backscattering and two extinction 
coefficients. The retrieval of particle microphysical parame- 
ters from these 3/6 + 2a data using inversion with regular- 
ization was discussed in our earlier publication, where good 
agreement between AERONET and lidar observations was 
reported (Veselovskii et al., 2009). 

Figure 4 shows the vertical profiles of aerosol backscat- 
tering and extinction coefficients measured at 355, 532 and 
1064nm wavelengths on 15 August 2006. The backscatter- 
ing shows a maximum at 1250 m and a secondary maximum 



Fig. 3. Uncertainty of particle volume estimation for different 
ranges of considered mp [0, 0.01], [0, 0.02]. [0, 0.03]. Simula- 
tion was performed for distribution with (a) ?'o = 0.2 pm and 
(b) i'q = 2 pm. Input errors are e = 10 % and model refractive index 
m = 1.45-10.005. 


at 1900 m. The particle size distribution for this day was rep- 
resented mainly by the fine mode and the uncertainties of 
the optical data measurements were estimated to be below 
10% (Veselovskii et al., 2009). The vertical profiles of vol- 
ume density, effective radius and real part of refractive in- 
dex obtained with the regularization and LE approaches are 
shown in Fig. 5 where the results obtained with both tech- 
niques are similar. The volume density profile has a sim- 
ilar shape as the particle backscattering, meaning that the 
particle radius and refractive index do not change signifi- 
cantly with height. The retrieved effective radius, as shown 
in Fig. 5b, is about 0.22 ± 0.055 pm for all heights, the uncer- 
tainty of retrieval is estimated from the results of numerical 
simulations summarized in Table 1 . It should be noted that 
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Fig. 4. Vertical profiles of aerosol backscattering (solid lines) and 
extinction (dashed lines) coefficients measured at 355, 532 and 1064 
nm wavelengths on 15 August 2006. 

the vertical profile of the effective radius obtained with LE 
oscillates less than the profile obtained with regularization, 
suggesting a more stable inversion. The refractive indices re- 
trieved with both techniques agree reasonably well. The real 
part of the refractive index slightly rises with height from 
1.37 ±0.05 to 1.43 ±0.05, the imaginary part m\ is below 
0.005 for all heights. 

As discussed in the previous section, the number of in- 
put optical data can be reduced when only bulk particle 
properties are desired. To test this claim, we also performed 
the inversion using the reduced set of optical data given by 
3/3 + la, where extinction at 532 nm is removed. The cor- 
responding results are also shown in Fig. 5. The inversion 
using either the full (3/3 + 2a) or reduced (3/1+ la) data 
sets leads to similar results, supporting the conclusions made 
from the numerical simulations. This comparison of the reg- 
ularization and LE approaches illustrates that the LE tech- 
nique can provide trustworthy estimations of particle param- 
eters. At the same time, the high speed of the retrieval us- 
ing linear estimates makes the method preferable for gen- 
erating bulk aerosol information from long-term series of 
lidar observations. 

5 Inversion of long-term series of multiwavelength lidar 
observations 

To test the retrieval of time-sequences of particle parame- 
ters we used data from the multiwavelength Raman lidar at 
TUBITAK Research Center located in the vicinity of Istan- 
bul, Turkey. The lidar is based on a frequency-tripled Quan- 
tel Brilliant B Nd:YAG laser with 10 Hz repetition rate. The 
pulse energies at k = 355, 532 and 1064 nm are 200, 250 and 




Fig. 5. Vertical profiles of (a) particle volume density, (b) effective 
radius, (c) real part of refractive index retrieved with LE approach 
from 3/3 + 2a and 3/0 + lev data and with regularization approach 
from 3/3 + 2 a data. 


300 mJ, respectively. The backscattered light is collected by 
a 40-cm aperture Newtonian telescope inclined so that the 
elevation angle is 40 degrees to the horizontal. The outputs 
of the detectors are recorded at 7.5 m range resolution using 
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Fig. 6. Particle (a) depolarization ratio and (b) extinction at 355 nm 
measured near Istanbul on 20 May 2010. Vertical resolution is 30 m 
for depolarization and 150 m for extinction. 


a Licel data acquisition system that incorporates both analog 
and photon-counting electronics. The system is able to moni- 
tor backscattering at 355, 532, 1064 nm, Raman nitrogen sig- 
nals at 387, 608 nm and Raman water vapor signal at 408 nm. 
A polarizing beamsplitter cube in the 355 nm channel allows 
simultaneous monitoring of co- and cross-polarized compo- 
nents of backscattered radiation. The particle depolarization 
ratio was calculated from the ratio of co- and cross polarized 
components of the particle backscattering coefficients. For 
the calibration of depolarization measurements the molecu- 
lar depolarization ratio in an aerosol-free region was used. 
In each file 3000 laser pulses were accumulated, thus the 
temporal resolution of the measurements is 5 min. 

The measurements were performed during May 2010 
when weak ash layers from the Eyjafjallajokull volcanic 
eruption periodically reached Turkey. The temporal evolu- 
tion of the particle depolarization ratio <5 p at 355 nm during 
the night of 20-21 May is shown in Fig. 6 . The highly depo- 
larizing volcanic layer with maximum particle depolarization 
ratio of approximately <5 p = 20 ± 5 % appears at 22:00 UTC 
and is observed for about a period of approximately four 
hours at 2-3 km heights. During the same time but for al- 
titudes below 2 km, the particle depolarization ratio did not 
exceed 5 %. The depolarization ratio in the ash layer is lower 
than the values of <5 p ~ 40 % that were observed over North- 
ern Europe (Ansmann et al., 2010), implying that the ash may 
have been mixed with more locally produced aerosols. The 



Fig. 7. Time-series of (a) particle volume density; (b) effective ra- 
dius; (c) real part of refractive index retrieved from measurements 
on 20 May 2010. 


analysis of the meteorological situation for this day and the 
optical properties of ash particles is given in Papayannis et 
al. (2011). The same figure shows the temporal variation of 
the aerosol extinction at 355 nm. The extinction is calculated 
from the Raman nitrogen signal (Ansmann et al., 1992). To 
decrease the uncertainty we reduced the height resolution up 
to 200 m and introduced a 3 point sliding average in the tem- 
poral domain. We estimate the uncertainty of the particle ex- 
tinction and backscattering 0 : 355 , /I 355 , /I 532 , Pwm calculation 
at the heights of interest (below 2.5 km) to be less than 10 %. 
Comparing Fig. 6 a and b we conclude that the optical depth 
of the ash layer is quite low (below 0.05). The Angstrom ex- 
ponent calculated from the extinction coefficients at 355 and 
532 nm was about 1 .8 in 0.8-2 km height range, meaning that 
the particles were relatively small. 


www.atmos-meas-tech.net/5/1135/2012/ 


Atmos. Meas. Tech., 5, 1135-1145, 2012 





1144 


I. Veselovskii et al.: Linear estimation of particle bulk parameters 


To retrieve the time-sequences of particle parameters, we 
used the 3/1 + 1 a data set, because the uncertainty of extinc- 
tion at 532 nm was too high for the chosen temporal resolu- 
tion. It also allowed us to test the ability of a reduced dataset 
to provide useful time series results. Figure 7a shows the 
time-height distribution of the particle volume density, which 
is similar to the time-height extinction distribution in Fig. 6b. 
The region of enhanced volume density is contoured and 
shown also in Fig. 6a in order to illustrate that it coincides 
well with the region of enhanced particle extinction. This im- 
plies that the particle size and refractive index did not vary 
significantly in this region. In the color maps in Fig. 7b, c 
showing effective radius and the real part of refractive index, 
the regions where the particle extinction is low are removed, 
because no reliable retrieval could be performed there. The 
particle effective radius is about 0.22 pm in 0.8-2 km height 
range and it does not vary significantly over the night. Some 
increase of effective radius is observed near the ash plume. 
The retrievals inside the ash layer should be taken with care 
because ash particles are of irregular shape and the retrieval 
is based on Mie kernels, assuming a spherical particle shape, 
that introduces significant uncertainties. In particular the real 
part of refractive index is significantly underestimated with 
this approach (Veselovskii et al., 2010). For the treatment 
of non-spherical particles, the kernels corresponding to ran- 
domly oriented spheroids (Dubovik et al., 2006) can be im- 
plemented as previously shown (Veselovskii et al., 2010), but 
that effort goes beyond the scope of present paper. 

The real part of refractive index shown in Fig. 7c varies in 
the range of 1.39-1.45. The marked region is characterized 
by a value in r ~ 1.4, indicating that the aerosol contains a 
significant amount of water. At low altitudes after midnight 
some enhancement of m r up to 1.45 is observed. As men- 
tioned, above 2 km the real part of the refractive index can 
be underestimated due to particle non-sphericity. The imagi- 
nary part of refractive index was estimated as 0.006 ± 0.003. 
The enhancement of m\ up to 0.01 was observed inside the 
ash layer, but again, for accurate quantification of m\ in this 
layer the spheroidal model should be used. Thus obtained 
results look reasonable and demonstrate that the use of lin- 
ear estimates makes possible the fast retrieval of time-height 
distributions of particle parameters from extended lidar ob- 
servations. The inversion of optical data to the aerosol pa- 
rameters shown in Fig. 7 took approximately 5 min using a 
standard laptop computer illustrating the potential of the LE 
technique for processing large volumes of the MW lidar data. 

6 Conclusions 

An algorithm for the linear estimation of aerosol bulk prop- 
erties such as particle volume and complex refractive in- 
dex from multiwavelength lidar measurements is presented. 
The particle concentration is estimated from a linear combi- 
nation of aerosol backscattering and extinction coefficients 


measured by multi-wavelength lidar while avoiding the re- 
trieval of the particle size distribution. This approach is 
shown to both increase the speed and stability of the inver- 
sion. The definition of the coefficients required for the lin- 
ear estimates is based on an expansion of the particle size 
distribution in terms of the measurement kernels. Once the 
coefficients for the linear estimates are established, the ap- 
proach allows very fast retrieval of aerosol bulk properties. 
In addition, the straightforward estimation of bulk proper- 
ties stabilizes the inversion making it more resistant to noise 
in the optical data: the retrieval does not fail even for input 
random uncertainties as large as 30 %. The uncertainties of 
the retrieval derived from numerical simulations are close to 
the values reported previously for the full inversion scheme 
that was used to derive the entire family of solutions using 
the regularization procedure (Veselovskii et al., 2002, 2004). 
The application of both techniques to the same lidar mea- 
surements did not reveal significant differences in the results 
of the two retrieval approaches. 

An important finding of this study is that it is feasible to 
reduce the number of input optical characteristics and still 
retrieve useful bulk aerosol properties. A comparison of in- 
versions using 3/1 + 2 a and 3/1 + la data demonstrates that 
excluding particle extinction at 532 nm does not significantly 
degrade the retrieval. At the same time, removing extinction 
at 355 nm enhances uncertainties of retrieval 

The high speed of the retrieval using linear estimates 
makes the method preferable for generating aerosol in- 
formation from long-term series of lidar observations. To 
demonstrate the efficiency of the method long-term series 
of aerosol physical properties derived from lidar observa- 
tions performed in Turkey in May 2010 were processed. As 
a result, the multi-wavelength lidar data, for the first time, 
were inverted into time-height distributions of particle pa- 
rameters. We should mention though that the algorithm stud- 
ied here should not be considered as a replacement for the 
full inversion (regularization) approach, because in many ap- 
plications 3/1 + 2 a data exist and the retrieval of PSD is crit- 
ical. However, if the data set is reduced, the current work 
demonstrates clearly that useful physical information may 
still be retrievable. 
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